home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / tests / testqp.sci < prev   
Text File  |  1999-09-16  |  7KB  |  240 lines
  1. function [ident,info]=testqp(n,irand,idef,imi,icota,nul,eps,modo)
  2. //Syntax : [ident,info]=testqp(n,irand,idef,imi,icota,nul,eps,modo)
  3. //
  4. //    test problems for quapro and linpro
  5. //
  6. //%inputs
  7. //   n        : dimension of the problem
  8. //   irand    : seed for starting the generation of random numbers.
  9. //              For irand < 0 the seed will not be recovered
  10. //   idef     : it indicates the type of quadratic problem. It may
  11. //              take the following values:
  12. //                  0  : for a linear problem
  13. //                  1  : for an indefinite quadratic problem
  14. //                  2  : for a positive semi-definite quadratic problem
  15. //   imi      : biggest number of equality or inequality constraints
  16. //   icota    : it indicates the type of bound constraints. It may
  17. //              take the values:
  18. //                > 2  : there are only security bounds
  19. //                  2  : there are security bounds for some variables
  20. //                  1  : there is at least one real bound for each
  21. //                       variable
  22. //   nul      : number of null columns of the Hessian
  23. //   eps      : maximum absolute differences for predicting active
  24. //              inequality constraints in macro postqp (standard
  25. //              value=%eps**0.75)
  26. //   modo     : it indicates the way the process begins:
  27. //                  1  : no initial point has been given. The
  28. //                       program computes a feasible initial point
  29. //                       that is a vertex of the region of feasible
  30. //                       points
  31. //                  2  : no initial point has been given. The program
  32. //                       computes a feasible initial point that is not
  33. //                       necesarily a vertex. This value of  MODO is
  34. //                       advisable when the quadratic form is positive
  35. //                       definite and there are a few constraints in
  36. //                       the problem or when there are large bounds
  37. //                       on the variables that are security bounds and
  38. //                       very probably  no actives at the solution
  39. //                  3  : an initial point is given
  40. //
  41. //%outputs
  42. //   ident    : seed for re-starting the generation of random numbers
  43. //              when there are several problems to be tested
  44. //   info     : flag showing return conditions:
  45. //                  0  : normal ending. A local minimum has been found
  46. //                > 0  : abnormal ending
  47. //!
  48.   gigant=1.e3; xci=10; xcs=100; r=gigant/n; t1=r/4;
  49.   imp=1; ident=irand; x0=0;
  50.   if irand > 0 then,..
  51.       rand('seed');..
  52.       rand(ones(1,irand)); ..
  53.   end;
  54. //
  55. // Generation of H
  56. //
  57.   if idef=0 then,..
  58.      h(n,n)=0;..
  59.   else,..
  60.      t=rand; ident=ident+1; s=nul*n/(2*gigant);
  61.      if t < s then,..
  62.         t=rand*nul; mnul=int(t); ident=ident+1;..
  63.      else,..
  64.         mnul=0;..
  65.      end;
  66.      h=rand(ones(n,n));
  67.      h=h+h';
  68.      ident=ident+n**2;
  69.      if mnul <> 0 then,..
  70.         nk=n/mnul;..
  71.      else,..
  72.         nk=0;..
  73.      end;
  74.      i1=nk; i10=nk;
  75.      for j=1:n,..
  76.         if j=i1 then h(:,j)=0; end,..
  77.         i1=i1+nk;..
  78.      end;
  79.      for j=1:n,..
  80.         if j=i10 then h(j,:)=0; end;..
  81.         i10=i10+nk;..
  82.      end;
  83.      if idef=2 then h=h'*h; end;..
  84.   end;
  85. //
  86. // Generation of p
  87. //
  88.   p=rand(n,1); ident=ident+n;
  89. //
  90. // constraints
  91. //
  92.   t=rand*imi; mi=int(t); ident=ident+1;
  93.   t=rand*imi; md=int(t); ident=ident+1;
  94.   mid=mi+md;
  95.   c=rand(n,mid); d=rand(mid,1); ident=ident+mid*(n+1);
  96. //
  97. // bounds
  98. //
  99.   if icota <= 2 then,..
  100.      for i=1:n,..
  101.         t=rand; ident=ident+1;..
  102.         if t <= 0.25d0 then,..
  103.            if icota = 1 then,..
  104.               ci(i)=xci; cs(i)=xcs;..
  105.            else,..
  106.               ci(i)=-gigant; cs(i)=gigant;..
  107.            end;..
  108.         else,..
  109.            if t <= 0.5d0 then,..
  110.               ci(i)=(t-0.375d0)*t1+xci; cs(i)=gigant;..
  111.            else,..
  112.               if t <= 0.75d0 then,..
  113.                  ci(i)=-gigant; cs(i)=(t-0.625d0)*t1+xcs;..
  114.               else,..
  115.                  ci(i)=(t-0.875d0)*t1+xci; cs(i)=ci(i)+rand*r+xcs;..
  116.                  ident=ident+1;..
  117.               end;..
  118.            end;..
  119.         end;..
  120.      end,..
  121.   end;
  122.   if icota > 2 then,..
  123.    cs=gigant*ones(n,1) ; ci=-cs ; ..
  124.   end;
  125. //
  126. //Point initial
  127. //
  128.   if modo=3 then,..
  129.      x1=0; h1=eye(n,n); p1(n)=0; modo1=2; imp1=0;..
  130.      [x0,f,lagr,info]=quapro(x1,h1,p1,c,d,ci,cs,mi,modo1,imp1);..
  131.   end;
  132. //
  133. //
  134.   write(%io(2), '  ************** call linpro/quapro ');
  135.   if idef=0 then,..
  136.      [x,f,lagr,info]=linpro(x0,p,c,d,ci,cs,mi,modo,imp);..
  137.   else,..
  138.      [x,f,lagr,info]=quapro(x0,h,p,c,d,ci,cs,mi,modo,imp);..
  139.   end;
  140. //
  141.   write(%io(2), ' *************** check');
  142.   if info = 0 then,..
  143.      [n1qp,n2qp,n3qp,n4qp,nact]=postqp(x,lagr,h,p,c,d,ci,cs,mi,eps);..
  144.   else,..
  145.      write(%io(2),'INCORRECT END!!!. INFO:');..
  146.      write(%io(2),info,'(10x,f6.0)');..
  147.   end
  148.  
  149.  
  150. function [n1qp,n2qp,n3qp,n4qp,nact]=postqp(x,lagr,h,p,c,d,ci,cs,mi,eps)
  151. //Syntax : <n1qp,n2qp,n3qp,n4qp,nact>=postqp(x,lagr,h,p,c,d,ci,cs,mi,eps)
  152. //
  153. //Controle de la programmation quadratique
  154. //
  155. //%parametres d'entree
  156. //   x,lagr : output of quapro
  157. //   h,p,c,d,ci,cs,mi : input of quapro
  158. //   eps      : maximum absolute differences for predicting active
  159. //              inequality constraints (standar value=%eps**0.75)
  160. //
  161. //output
  162. //    n1qp : kuhn-tucker
  163. //    n2qp : error in constraints
  164. //    n3qp : lagrange
  165. //    n4qp : 
  166. //    nact : active constraints
  167. //
  168. //    n1qp n2qp n3qp and n4qp must be small
  169. //!
  170.   n=maxi(size(p));
  171.   slagr=size(lagr); md=slagr(1)-n-mi; nact=mi;
  172. //
  173. //n1qp=
  174.   mid=mi+md;
  175. //
  176.   n1qp=p + h*x + lagr(1:n);
  177.   if mid > 0 then n1qp=n1qp + c*lagr(n+1:n+mid);end;
  178.   n1qp=norm(n1qp)
  179. //
  180. //n2qp=
  181. //
  182.   n2qp=0;
  183.   if mid > 0 then,..
  184.      cxmd=c'*x-d;..
  185.   end;
  186.   if mi>0 then n2qp=norm(cxmd(1:mi),1);end;
  187.   for i=mi+1:mi+md,..
  188.      n2qp=n2qp+maxi(cxmd(i),0);..
  189.   end;
  190. //
  191.   for i=1:n,..
  192.      n2qp=n2qp+maxi(ci(i)-x(i),0)+maxi(0,x(i)-cs(i));..
  193.   end;
  194. //
  195. //  n3qp=
  196. //  n4qp=
  197. //
  198.   n3qp=maxi(abs(lagr));
  199.   n4qp=0;
  200.   for i=1:n,..
  201.      s1=abs(x(i)-ci(i));..
  202.      if s1 < eps then,..
  203.         n3qp=mini(n3qp,-lagr(i));..
  204.         n4qp=n4qp+abs(lagr(i)*(x(i)-ci(i)));..
  205.         nact=nact+1;..
  206.      end,..
  207.      s2=abs(x(i)-cs(i));..
  208.      if s2 < eps then,..
  209.         n3qp=mini(n3qp,lagr(i));..
  210.         n4qp=n4qp+abs(lagr(i)*(x(i)-cs(i)));..
  211.         nact=nact+1;..
  212.      end,..
  213.   end;
  214. //
  215.   if md > 0  then,..
  216.      for i=1:md,..
  217.         if abs(cxmd(mi+i)) < eps then,..
  218.            n3qp=mini(n3qp,lagr(n+mi+i));..
  219.            nact=nact+1;..
  220.         end,..
  221.      end;..
  222.   end;
  223.   if mid > 0 then,..
  224.      n4qp=n4qp+abs(lagr(1+n:n+mid)'*cxmd(1:mid));..
  225.   end;
  226.   write(%io(2),'-NUMBER ACTIVE CONSTRAINTS');
  227.   write(%io(2),nact,'(10x,f3.0)');
  228.   write(%io(2),'-NORM OF  KUHN-TUCKER VECTOR:');
  229.   write(%io(2),n1qp,'(10x,e14.8)');
  230.   write(%io(2),'-ERROR IN CONSTRAINTS:');
  231.   write(%io(2),n2qp,'(10x,e14.8)');
  232.   write(%io(2),'-SMALLEST LAGRANGE MULTIPLIER:');
  233.   write(%io(2),n3qp,'(10x,e14.8)');
  234.   write(%io(2),'-COMPLEMENTARITY: multiplier*(c*x-d):');
  235.   write(%io(2),n4qp,'(10x,e14.8)')
  236.  
  237.  
  238.  
  239.  
  240.